REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway, 
Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection 
of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

08-31-2013  Final 

4.  TITLE  AND  SUBTITLE 


3.  DATES  COVERED  (From  —  To) 
1  June  2010  —  31  May  2013 

5a.  CONTRACT  NUMBER 


(YIP-10)  EFFICIENT  AND  ROBUST  HIGH-ORDER  METHODS  FOR 
FLUID  AND  SOLID  MECHANICS 


5b.  GRANT  NUMBER 


FA9550-10-1-0229 

5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


5 cl.  PROJECT  NUMBER 


Persson,  Per-Olof 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Department  of  Mathematics 
University  of  California,  Berkeley 
Berkeley,  CA  94720-3840 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFOSR/RSL 

875  North  Randolph  Street 
Suite  325,  Room  3112 
Arlington,  VA  22203 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approval  for  public  release;  distribution  is  unlimited. 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFOSR 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

The  goal  of  the  project  was  to  develop  new  numerical  schemes  and  solvers  for  high-order  accurate  simulations  of 
problems  in  fluid  and  solid  mechanics.  Three  main  areas  were  addressed  -  the  high  computational  cost  of  existing 
methods,  their  lack  of  robustness,  and  the  need  for  new  high-order  mesh  manipulation  algorithms.  The  project  has  led  to 
significant  developments  in  so-called  Line-DG  and  IMEX  based  numerical  schemes,  new  formulations  for  problems  with 
deforming  domains,  artificial  viscosity  based  stabilization,  and  mesh  deformation  using  a  nonlinear  elasticity  analogy. 

The  results  were  applied  to  a  number  of  important  real-world  problems,  such  as  drag  prediction  for  turbulent  flows, 
flapping  flight  with  fluid-structure  interaction,  aeroacoustics  problems,  and  transonic/supersonic  flow  problems.  The 
findings  were  disseminated  through  a  wide  range  of  publications,  presentations,  and  public  domain  software. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 


b.  ABSTRACT!  c.  THIS  PAGE 


17.  LIMITATION  OF  18.  NUMBER  19a.  NAME  OF  RESPONSIBLE  PERSON 
ABSTRACT  PAGES  ^T'  Per"01of  Persson 

19b.  TELEPHONE  NUMBER  (include  area  code ) 

UU  1  (510)  642-6947 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


FINAL  REPORT 


(YIP-10)  EFFICIENT  AND  ROBUST  HIGH-ORDER  METHODS 
FOR  FLUID  AND  SOLID  MECHANICS 

Per-Olof  Persson 

Department  of  Mathematics 
University  of  California,  Berkeley 
Berkeley,  CA  94720-3840 

August  2013 


Grant  number:  FA9550-10-1-0229 
Dates  covered:  1  June  2010  —  31  May  2013 


Abstract 

The  goal  of  the  project  was  to  develop  new  numerical  schemes  and  solvers  for  high-order  accu¬ 
rate  simulations  of  problems  in  fluid  and  solid  mechanics.  Three  main  areas  were  addressed  - 
the  high  computational  cost  of  existing  methods,  their  lack  of  robustness,  and  the  need  for  new 
high-order  mesh  manipulation  algorithms.  The  project  has  led  to  significant  developments  in 
so-called  Line-DG  and  IMEX  based  numerical  schemes,  new  formulations  for  problems  with  de¬ 
forming  domains,  artificial  viscosity  based  stabilization,  and  mesh  deformation  using  a  nonlinear 
elasticity  analogy.  The  results  were  applied  to  a  number  of  important  real-world  problems,  such 
as  drag  prediction  for  turbulent  flows,  flapping  flight  with  fluid- structure  interaction,  aeroacous- 
tics  problems,  and  transonic/supersonic  flow  problems.  The  findings  were  disseminated  through 
a  wide  range  of  publications,  presentations,  and  public  domain  software. 


Efficient  and  Robust  High-Order  Methods  for  Fluid  and  Solid  Mechanics 

P.-O.  Persson,  Dept,  of  Mathematics,  University  of  California,  Berkeley 

Objectives 

The  goal  of  the  project  was  to  develop  new  numerical  schemes  and  solvers  for  high-order  accu¬ 
rate  simulations  of  problems  in  fluid  and  solid  mechanics.  Three  main  areas  were  addressed  -  the 
high  computational  cost  of  existing  methods,  their  lack  of  robustness,  and  the  need  for  new  high- 
order  mesh  manipulation  algorithms.  The  project  has  led  to  significant  developments  in  so-called 
Line-DG  and  IMEX  based  numerical  schemes,  new  formulations  for  problems  with  deforming  do¬ 
mains,  artificial  viscosity  based  stabilization,  and  mesh  deformation  using  a  nonlinear  elasticity 
analogy.  The  results  were  applied  to  a  number  of  important  real-world  problems,  such  as  drag  pre¬ 
diction  for  turbulent  flows,  flapping  flight  with  fluid-structure  interaction,  aeroacoustics  problems, 
and  transonic/supersonic  flow  problems.  The  findings  were  disseminated  through  a  wide  range  of 
publications,  presentations,  and  public  domain  software. 


Accomplishments 

During  the  award  period,  we  have  conducted  several  significant  research  efforts  in  the  general  area 
of  high-order  accurate  numerical  methods  for  solving  partial  differential  equations  on  unstructured 
meshes.  It  is  widely  believed  that  these  new  methods  will  eventually  replace  the  more  traditional 
simulation  techniques  for  challenging  problems  in  fluid  and  solid  mechanics,  for  example  with 
propagating  waves,  turbulent  flow,  nonlinear  interactions,  and  multiple  scales.  However,  several 
new  developments  are  required  to  make  high-order  methods  competitive,  such  as  more  efficient 
numerical  solvers  and  increased  robustness.  Our  work  has  addressed  many  of  these  issues,  and 
below  we  summarize  the  research  that  we  have  carried  out. 


The  Line  Discontinuous  Galerkin  (Line-DG)  Method 


In  [8,  11],  we  presented  a  new  line-based  discontinuous  Galerkin  (DG)  discretization  scheme 
for  first-  and  second-order  systems  of  partial  differential  equations.  The  scheme  is  based  on  fully 
unstructured  meshes  of  quadrilateral  or  hexahedral  elements,  and  it  is  closely  related  to  the  standard 
nodal  DG  scheme  as  well  as  several  of  its  variants  such  as  the  collocation-based  DG  spectral  element 
method  (DGSEM)  or  the  spectral  difference  (SD)  method.  However,  our  motivation  is  to  maximize 
the  sparsity  of  the  Jacobian  matrices,  since  this  directly  translates  into  higher  performance  in 
particular  for  implicit  solvers,  while  maintaining  many  of  the  good  properties  of  the  DG  scheme. 

The  discretization  starts  by  mapping  a  first-order  system  of  conservation  laws  to  a  reference 
frame  in  terms  of  the  contravariant  fluxes,  see  Figure  1.  Each  of  the  spatial  derivatives  is  then 
discretized  separately  using  a  standard  DG  scheme: 
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where  the  numerical  contravariant  fluxes  can  be  written  in  terms  of  a  standard  approximate  Rie- 
mann  solver  for  the  original  fluxes,  since 

/i  =  F  ■  JV+  =  (JG~lF)  -N+  =  F  ( JG~TN+ )  =  F  ■  n+  (2) 
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Figure  1:  Top:  A  two-dimensional  illustration  of  the  mapping  from  a  reference  element  V  to  the 
actual  curved  element  v,  for  the  case  p  =  4.  Bottom:  The  number  of  connectivities  per  node  for 
a  first-order  operator  and  3-D  hexahedral  elements  with  the  Line-DG,  the  DGSEM/SD,  and  the 
nodal  DG  methods. 

Note  that  each  derivative  is  essentially  a  one-dimensional  discretization,  which  results  in  a  very 
high  level  of  sparsity.  The  resulting  scheme  is  similar  to  a  collocation  scheme,  but  it  uses  fully 
consistent  integration  along  each  1-D  coordinate  direction  which  results  in  different  properties  for 
nonlinear  problems  and  curved  elements.  Also,  the  scheme  uses  solution  points  along  each  element 
face,  which  further  reduces  the  number  of  connections  with  the  neighboring  elements.  This  sparsity 
pattern  is  illustrated  in  Figure  1  where  we  note  that  in  three  dimensions,  already  for  p  =  3  the  Line- 
DG  method  is  5.5  times  sparser  than  nodal  DG,  and  for  p  =  10  it  is  almost  40  times  sparser.  This 
reduction  in  stencil  size  translates  into  lower  assembly  times,  but  more  importantly,  for  matrix- 
based  solvers  it  means  drastically  lower  storage  requirements  and  faster  matrix- vector  products  for 
iterative  implicit  solvers.  For  second-order  derivatives,  we  use  an  LDG-type  scheme  with  consistent 
switches  along  globally  connected  lines. 

The  scheme  appears  well  suited  for  a  wide  range  of  problems,  and  we  have  demonstrated  optimal 
convergence  for  viscous  and  convective  model  problems,  as  well  as  the  full  compressible  Navier- 
Stokes  equations.  An  example  of  a  turbulent  flow  is  shown  in  figure  2,  where  the  Line-DG  method 
with  ILES  modeling  is  used  to  predict  the  acoustic  response  of  a  recorder. 

Efficient  LES  Time- Integration  using  Implicit- Explicit  Runge-Kutta  Schemes 

For  many  flow  problems  modeled  by  Large  Eddy  Simulation  (LES),  the  computational  meshes 
are  such  that  a  majority  of  the  elements  would  allow  for  explicit  timestepping,  but  the  CFL- 
condition  is  limited  by  the  smaller  stretched  elements  in  the  boundary  layers.  We  have  developed 
an  implicit-explicit  time-integration  scheme  that  uses  an  implicit  solver  only  for  the  smaller  portion 
of  the  domain  where  it  is  required  to  avoid  severe  timestep  restrictions,  but  an  efficient  explicit  solver 
for  the  rest  of  the  domain  [2].  We  use  the  Runge-Kutta  IMEX  schemes  and  have  demonstrated  the 
stability  and  accuracy  of  the  solver  for  a  wide  range  of  schemes  and  orders  of  accuracy.  In  addition, 
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Figure  2:  Simulation  of  the  aeroacoustics  of  a  recorder,  using  the  Line-DG  scheme  and  ILES 
modeling.  The  left  plots  show  the  mesh,  with  polynomial  degrees  p  =  7  within  each  element,  and 
the  right  plot  shows  a  sample  solution  (vorticity). 
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Figure  3:  Splitting  of  the  mesh  around  an  SD7003  airfoil  into  implicit  and  explicit  elements.  Less 
than  9%  of  the  elements  are  in  the  boundary  layer  and  integrated  using  the  implicit  scheme.  The 
right  plot  shows  an  instantaneous  flow  solution,  as  isosurfaces  of  the  entropy  colored  by  the  Mach 
number. 

we  use  a  quasi-Newton  method  for  the  implicit  region,  which  allows  us  to  re-use  the  Jacobian 
matrices  and  their  incomplete  factorizations  for  a  large  number  of  timesteps.  We  also  showed  the 
application  of  the  technique  to  a  realistic  LES-type  problem  of  turbulent  flow  around  an  airfoil, 
where  we  conclude  that  the  approach  can  give  performance  that  is  superior  to  both  fully  explicit 
and  fully  implicit  methods  (see  figure  3).  The  histogram  of  the  elements  sizes  shows  how  about  9% 
of  the  elements  are  treated  implicitly. 

IMEX  schemes  for  Fluid- Structure  Interaction 

Many  important  scientific  and  engineering  problems  require  predictions  of  fluid-structure  inter¬ 
action  (FSI).  For  example,  oscillatory  interactions  in  engineering  systems  (e.g.  aircraft,  turbines, 
and  bridges)  can  lead  to  failure.  The  blood  flow  in  arteries  and  artificial  heart  valves  is  highly 
dependent  on  structural  interactions.  These  interactions  often  involve  multiple  scales  and  non¬ 
linear  effects,  which  makes  it  challenging  to  solve  even  relatively  simple  problems  accurately.  In 
[17,  19],  we  presented  a  high-order  accurate  scheme  for  fully  coupled  FSI  problems.  Using  implicit- 
explicit  (IMEX)  Runge-Kutta  methods,  we  treat  a  predicted  traction  t  explicitly  and  everything 
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Figure  4:  Fluid-structure  interaction  simulation  of  the  flow  around  a  three  dimensional  flexible 
membrane  at  various  times  (Mach  number  on  iso-entropy  surfaces).  Compared  to  a  rigid  plate,  the 
leading  edge  of  the  membrane  aligns  with  the  fluid  and  reduces  the  leading  edge  separation.  The 
right  plot  shows  a  sample  deformed  mesh  used  in  the  ALE  formulation  for  the  fluid  problem. 


else  implicitly: 
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where  the  superscripts  /  and  s  correspond  to  the  fluid  and  the  structure,  respectively.  This  leads 
to  a  block  upper  triangular  system  that  fully  decouples  the  implicit  solution  procedures  for  the 
fluid  and  the  solid  parts,  which  we  perform  using  two  separate  efficient  parallel  solvers.  We  also 
demonstrated  using  a  model  problem  that  subiterations  are  generally  not  required,  which  gives  the 
scheme  close  to  the  optimal  performance  of  solving  the  fluid  and  structure  separately.  An  example 
of  highly  separated  flow  around  a  thin  membrane  is  shown  in  Figure  4. 

Artificial  viscosity  based  shock  capturing  for  transient  flow  problems 

We  have  previously  introduced  a  widely  popular  artificial  viscosity  based  scheme  for  sub-cell 
shock  capturing  with  high-order  DG  schemes.  A  highly  selective  sensor  is  based  on  the  decay  rate  of 
the  expansion  coefficients  of  the  solution  in  an  orthogonal  basis.  The  scheme  is  an  attractive  alter¬ 
native  to  limiting,  partly  because  of  the  subgrid  accuracy,  the  ability  to  obtain  Newton-converged 
solutions,  and  the  ease  of  implementation.  We  have  recently  continued  this  work  and  developed 
similar  schemes  for  transient  problems  with  moving  shocks,  time-stepped  using  high-order  accurate 
implicit  schemes  with  Jacobian-based  Newton-Krylov  solvers.  In  [16],  we  demonstrated  that  the 
sensors  can  be  coupled  weakly  without  losing  robustness,  which  simplifies  the  implementation  and 
reduces  the  computational  times.  The  weak  coupling  also  allows  for  non-compact  regularization 
of  the  artificial  viscosity,  and  we  show  that  C°  and  C 2  regularity  in  the  sensor  greatly  improves 
the  solution  quality.  The  hypersonic  double  mach  reflection  problem  of  Woodward  and  Colella  is 
shown  in  figure  5  (top),  as  the  density  field  and  the  applied  artificial  viscosity.  In  the  bottom  plot 
we  show  a  3-D  problem  with  the  transonic  flow  around  a  tapered  wing,  as  pressure  and  the  applied 
artificial  viscosity.  In  both  cases,  the  plots  clearly  indicate  how  the  sensor  is  applied  only  to  a 
narrow  band  of  elements. 

Arbitrary  Lagrangian-Eulerian  Formulations 

Several  approaches  have  been  proposed  for  handling  CFD  simulations  that  involve  time- varying 
domains  with  moving  or  deformable  bodies.  One  particularly  successful  method  is  the  Arbitrary 
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Figure  5:  Top:  The  double-mach  reflection  problem  of  Woodward  and  Colella,  showing  density 
(left)  and  the  C°-continuous  artificial  viscosity  (right).  Bottom:  The  transonic  flow  around  a 
tapered  wing,  pressure  (left)  and  the  artificial  viscosity  (right). 

Lagrangian-Eulerian  (ALE)  formulation,  where  the  mesh  is  deformed  arbitrarily  but  conforming  to 
the  moving  boundaries.  The  physical  equations  are  modified  with  convective  terms  to  account  for 
the  grid  motion.  We  have  developed  a  mapping-based  ALE  scheme  where  the  modified  equations 
are  always  solved  on  a  fixed  grid.  This  has  the  advantage  that  the  scheme  is  stable  and  accurate  for 
arbitrary  orders  of  accuracy  in  space  and  time.  The  so-called  Geometric  Conservation  Law  (GCL) 
is  handled  by  solving  an  additional  equations  for  the  Jacobian  of  the  mapping.  Using  a  nonlinear 
elasticity  approach,  we  can  produce  mappings  for  highly  stretched  and  deformed  domains,  see 
Figure  6  for  an  example  of  a  flapping  bat  and  a  corresponding  flow  field.  We  have  also  derived 
so-called  stage-consistent  mesh  velocities  for  implicit  Runge-Kutta  schemes  [19],  which  is  a  general 
framework  for  deriving  high-order  accurate  numerical  mappings  from  a  set  of  deformed  meshes. 

High  Fidelity  Simulation  of  Optimized  Flapping  Wing  MAVs 

In  [3,  7,  14],  we  proposed  a  multi-fidelity  framework  for  the  design  of  flapping  wing  Micro  Aerial 
Vehicles.  Since  high-fidelity  solvers  are  prohibitively  expensive  for  using  in  the  optimization  process, 
we  use  a  wake-only  method  and  a  fast  panel  code  to  find  candidate  designs,  which  are  validated 
using  our  high-order  Navier-Stokes  code.  The  deforming  domain  is  discretized  using  our  mapping- 
based  ALE  approach.  The  need  for  high-order  methods  is  clearly  seen  in  our  p-convergence  studies, 
where  the  forces  vary  substantially  with  polynomial  orders  p  —  1  and  p  —  2  but  essentially  show 
grid  convergence  for  p  —  3  and  higher. 

In  our  design  study,  we  found  that  the  designs  obtained  by  the  low-fidelity  solvers  often  perform 
poorly  when  including  the  viscous  effects  of  the  high-fidelity  DG  solver,  due  to  the  large  amounts  of 
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Figure  6:  Our  high-order  mapping-based  ALE  approach  [7,  19]  for  the  simulation  of  the  flapping 
flight  of  a  bat.  The  reference  mesh  (left)  is  deformed  in  time  using  a  non-linear  elasticity  approach 
which  maintains  the  high-quality  of  the  elements  (middle).  The  right  plot  shows  a  sample  solution 
field. 


Figure  7:  Four  snapshots  in  time  of  the  optimal  flapping  wing  simulation,  without  camber  (top)  and 
with  camber  (bottom) .  The  flow  separation  is  significantly  reduced  by  using  a  dynamic  cambering. 


separation  resulting  from  the  wing  twisting  (Figure  7,  top  row).  To  reduce  this  effect,  we  introduce 
a  dynamically  varying  camber  along  the  span,  which  is  aligning  the  leading  edge  with  the  incoming 
flow.  This  result  in  a  flow  that  stays  attached  for  most  of  the  flapping  cycle,  except  for  some  minor 
tip  vortices,  and  ultimately  a  higher  performance  of  the  flapping  wing  pair  (Figure  7,  bottom  row). 
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